f = oda ~ oda_lag + p_polity2_lag + gdppc_lag + pop_lag + factor(year)*factor(agency)
f1 = update(f, '.~.+alliance_lag*independence_lag')
f2 = update(f, '.~.+trade_lag*independence_lag')
f3 = update(f, '.~.+infant_mortality_lag*independence_lag')

mod1 = lm(f1, dat)
mod2 = lm(f2, dat)
mod3 = lm(f3, dat)

lags = list(mod1, mod2, mod3)
se = lapply(lags, function(i) sqrt(diag(sandwich(i))))

print(screenreg(lags, omit.coef='factor', digits=3, override.se=se))
